################################################
# Analysis for:
#
# Warner, Zach, Garrett N. Vande Kamp, and Soren Jordan.
# ``Conditional Relationships in Dynamic Models.'' 
#  Political Science Research & Methods
#
# Cavari Supplemental Analysis (SM Section 6; Section 8.1)
#
# Required files: main-cavari-psq
#                 interp_urdf.R (for interpretation)
#
################################################

# Set your working directory here

library(dplyr)
library(ggplot2)
library(grid)
library(haven)
library(lmtest)
library(dynamac)
library(urca)
library(tseries)
library(polywog)
library(MASS)

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  if (is.null(layout)) {
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    for (i in 1:numPlots) {
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

# ur.df output is annoying. interp_urdf.R helper function for interpretation
#  Author: Hank Roark. 
#  https://gist.github.com/hankroark/968fc28b767f1e43b5a33b151b771bf9
source("interp_urdf.R")

##############################################################################################################################
# In-text replication: Cavari (2019)
##############################################################################################################################
approval <- read_dta("main-cavari-psq.dta")

# Estimating the Original Model
approval$rr <- ifelse(approval$PRESIDENT == 40, 1, 0)
approval$ghwb <- ifelse(approval$PRESIDENT == 41, 1, 0)
approval$bc <- ifelse(approval$PRESIDENT == 42, 1, 0)
approval$gwb <- ifelse(approval$PRESIDENT == 43, 1, 0)
approval$bo <- ifelse(approval$PRESIDENT == 44, 1, 0)

approval$ECONAPP_ECONMIP <- approval$APPROVE_ECONOMY*approval$MIP_MACROECONOMICS
approval$FPAPP_FPMIP <- approval$APPROVE_FOREIGN*approval$MIP_FOREIGN

# Create matrix of x regressors
approval.small <- approval %>%
				dplyr::select(APPROVE, APPROVE_ECONOMY, APPROVE_FOREIGN, MIP_MACROECONOMICS, MIP_FOREIGN, ECONAPP_ECONMIP, FPAPP_FPMIP, 
					APPROVE_L1, PARTY_IN, PARTY_OUT, UNRATE, DIVIDEDGOV, ELECTION, HONEYMOON, ghwb, bc, gwb, bo)

cavari.model <- lm(as.formula(paste0("APPROVE ~ ", paste(names(approval.small), collapse = "+"))), data = approval)

######### Table SM. 2; ``Cavari'' model #########
summary(cavari.model)
bgtest(cavari.model, order = 5, type = "F")


## Extend to the general model

# Creating Lags 
approval$MIP_FOREIGN_L1 <- lshift(approval$MIP_FOREIGN, 1)
approval$MIP_MACROECONOMICS_L1 <- lshift(approval$MIP_MACROECONOMICS, 1)
approval$APPROVE_ECONOMY_L1 <- lshift(approval$APPROVE_ECONOMY, 1)
approval$APPROVE_FOREIGN_L1 <- lshift(approval$APPROVE_FOREIGN, 1)
approval$UNRATE_L1 <- lshift(approval$UNRATE, 1)
approval$PARTY_IN_L1 <- lshift(approval$PARTY_IN, 1)
approval$PARTY_OUT_L1 <- lshift(approval$PARTY_OUT, 1)

approval$ECONAPP_L1_ECONMIP_L1 <- approval$APPROVE_ECONOMY_L1*approval$MIP_MACROECONOMICS_L1
approval$FPAPP_L1_FPMIP_L1 <- approval$APPROVE_FOREIGN_L1*approval$MIP_FOREIGN_L1

approval$ECONAPP_ECONMIP_L1 <- approval$APPROVE_ECONOMY*approval$MIP_MACROECONOMICS_L1
approval$FPAPP_FPMIP_L1 <- approval$APPROVE_FOREIGN*approval$MIP_FOREIGN_L1

approval$ECONAPP_L1_ECONMIP <- approval$APPROVE_ECONOMY_L1*approval$MIP_MACROECONOMICS
approval$FPAPP_L1_FPMIP <- approval$APPROVE_FOREIGN_L1*approval$MIP_FOREIGN

#Testing stationarity of all variables and interactions
df_ECONAPP_ECONMIP <- ur.df(approval$ECONAPP_ECONMIP, type = "drift", selectlags = "AIC")
interp_urdf(df_ECONAPP_ECONMIP)
df_FPAPP_FPMIP <- ur.df(approval$FPAPP_FPMIP, type = "drift", selectlags = "AIC")
interp_urdf(df_FPAPP_FPMIP)

df_ECONAPP_ECONMIP_L1 <- ur.df(approval$ECONAPP_ECONMIP_L1[2:140], type = "drift", selectlags = "AIC")
interp_urdf(df_ECONAPP_ECONMIP_L1)
df_FPAPP_FPMIP_L1 <- ur.df(approval$FPAPP_FPMIP_L1[2:140], type = "drift", selectlags = "AIC")
interp_urdf(df_FPAPP_FPMIP_L1)

df_ECONAPP_L1_ECONMIP <- ur.df(approval$ECONAPP_L1_ECONMIP[2:140], type = "drift", selectlags = "AIC")
interp_urdf(df_ECONAPP_ECONMIP)
df_FPAPP_L1_FPMIP <- ur.df(approval$FPAPP_L1_FPMIP[2:140], type = "drift", selectlags = "AIC")
interp_urdf(df_FPAPP_L1_FPMIP)

df_APPROVAL <- ur.df(approval$APPROVE, type = "drift", selectlags = "AIC")
interp_urdf(df_APPROVAL)

df_APPROVE_ECONOMY <- ur.df(approval$APPROVE_ECONOMY, type = "drift", selectlags = "AIC")
interp_urdf(df_APPROVE_ECONOMY)
df_APPROVE_FOREIGN <- ur.df(approval$APPROVE_FOREIGN, type = "drift", selectlags = "AIC")
interp_urdf(df_APPROVE_FOREIGN)
df_MIP_MACROECONOMICS <- ur.df(approval$MIP_MACROECONOMICS, type = "drift", selectlags = "AIC")
interp_urdf(df_MIP_MACROECONOMICS) #Unit root!!!
df_MIP_FOREIGN <- ur.df(approval$MIP_FOREIGN, type = "drift", selectlags = "AIC")
interp_urdf(df_MIP_FOREIGN) #Unit root!!!
df_PARTY_IN <- ur.df(approval$PARTY_IN[1:138], type = "drift", selectlags = "AIC")
interp_urdf(df_PARTY_IN) #Unit root!!!
df_PARTY_OUT <- ur.df(approval$PARTY_OUT[1:138], type = "drift", selectlags = "AIC")
interp_urdf(df_PARTY_OUT)
df_UNRATE <- ur.df(approval$UNRATE, type = "drift", selectlags = "AIC")
interp_urdf(df_UNRATE) #Unit root!!!
df_DIVIDEDGOV <- ur.df(approval$DIVIDEDGOV, type = "drift", selectlags = "AIC")
interp_urdf(df_DIVIDEDGOV)

###Difference Unit Roots
approval$UNRATE_D <- approval$UNRATE - approval$UNRATE_L1
approval$MIP_MACROECONOMICS_D <- approval$MIP_MACROECONOMICS - approval$MIP_MACROECONOMICS_L1
approval$MIP_FOREIGN_D <- approval$MIP_FOREIGN - approval$MIP_FOREIGN_L1

# Create lags
approval$MIP_FOREIGN_D_L1 <- lshift(approval$MIP_FOREIGN_D, 1)
approval$MIP_MACROECONOMICS_D_L1 <- lshift(approval$MIP_MACROECONOMICS_D, 1)
approval$UNRATE_D_L1 <- lshift(approval$UNRATE_D, 1)

# Create lagged and cross-lagged interactions for differenced variable
# t/t
approval$ECONAPP_ECONMIP_D <- approval$APPROVE_ECONOMY*approval$MIP_MACROECONOMICS_D
approval$FPAPP_FPMIP_D <- approval$APPROVE_FOREIGN*approval$MIP_FOREIGN_D

# t-1/t-1
approval$ECONAPP_L1_ECONMIP_D_L1 <- approval$APPROVE_ECONOMY_L1*approval$MIP_MACROECONOMICS_D_L1
approval$FPAPP_L1_FPMIP_D_L1 <- approval$APPROVE_FOREIGN_L1*approval$MIP_FOREIGN_D_L1

# t/t-1
approval$ECONAPP_ECONMIP_D_L1 <- approval$APPROVE_ECONOMY*approval$MIP_MACROECONOMICS_D_L1
approval$FPAPP_FPMIP_D_L1 <- approval$APPROVE_FOREIGN*approval$MIP_FOREIGN_D_L1

# t-1/t
approval$ECONAPP_L1_ECONMIP_D <- approval$APPROVE_ECONOMY_L1*approval$MIP_MACROECONOMICS_D
approval$FPAPP_L1_FPMIP_D <- approval$APPROVE_FOREIGN_L1*approval$MIP_FOREIGN_D

approval_standard <- approval %>%
						filter(!is.na(UNRATE_D_L1), !is.na(PARTY_IN))


######### Table SM. 6 #########
adf.test(approval_standard$APPROVE_ECONOMY)$p.value
adf.test(approval_standard$APPROVE_FOREIGN)$p.value
adf.test(approval_standard$MIP_MACROECONOMICS)$p.value
adf.test(approval_standard$MIP_FOREIGN)$p.value
adf.test(approval_standard$MIP_MACROECONOMICS_D)$p.value
adf.test(approval_standard$MIP_FOREIGN_D)$p.value
adf.test(approval_standard$ECONAPP_ECONMIP_D)$p.value
adf.test(approval_standard$FPAPP_FPMIP_D)$p.value
adf.test(approval_standard$PARTY_IN)$p.value
adf.test(approval_standard$PARTY_OUT)$p.value
adf.test(approval_standard$UNRATE)$p.value
adf.test(approval_standard$UNRATE_D)$p.value
adf.test(approval_standard$DIVIDEDGOV)$p.value
adf.test(approval_standard$ELECTION)$p.value
adf.test(approval_standard$HONEYMOON)$p.value


Box.test(approval_standard$APPROVE_ECONOMY)$p.value
Box.test(approval_standard$APPROVE_FOREIGN)$p.value
Box.test(approval_standard$MIP_MACROECONOMICS)$p.value
Box.test(approval_standard$MIP_FOREIGN)$p.value
Box.test(approval_standard$MIP_MACROECONOMICS_D)$p.value
Box.test(approval_standard$MIP_FOREIGN_D)$p.value
Box.test(approval_standard$ECONAPP_ECONMIP_D)$p.value
Box.test(approval_standard$FPAPP_FPMIP_D)$p.value
Box.test(approval_standard$PARTY_IN)$p.value
Box.test(approval_standard$PARTY_OUT)$p.value
Box.test(approval_standard$UNRATE)$p.value
Box.test(approval_standard$UNRATE_D)$p.value
Box.test(approval_standard$DIVIDEDGOV)$p.value
Box.test(approval_standard$ELECTION)$p.value
Box.test(approval_standard$HONEYMOON)$p.value


kpss.test(approval_standard$APPROVE_ECONOMY)$p.value
kpss.test(approval_standard$APPROVE_FOREIGN)$p.value
kpss.test(approval_standard$MIP_MACROECONOMICS)$p.value
kpss.test(approval_standard$MIP_FOREIGN)$p.value
kpss.test(approval_standard$MIP_MACROECONOMICS_D)$p.value
kpss.test(approval_standard$MIP_FOREIGN_D)$p.value
kpss.test(approval_standard$ECONAPP_ECONMIP_D)$p.value
kpss.test(approval_standard$FPAPP_FPMIP_D)$p.value
kpss.test(approval_standard$PARTY_IN)$p.value
kpss.test(approval_standard$PARTY_OUT)$p.value
kpss.test(approval_standard$UNRATE)$p.value
kpss.test(approval_standard$UNRATE_D)$p.value
kpss.test(approval_standard$DIVIDEDGOV)$p.value
kpss.test(approval_standard$ELECTION)$p.value
kpss.test(approval_standard$HONEYMOON)$p.value

######### Table SM. 2; ``General'' model #########
cavari.model.differenced <- lm(APPROVE ~ APPROVE_ECONOMY + APPROVE_FOREIGN + MIP_MACROECONOMICS_D + MIP_FOREIGN_D +
                         ECONAPP_ECONMIP_D + FPAPP_FPMIP_D + 
                         APPROVE_ECONOMY_L1 + APPROVE_FOREIGN_L1 + MIP_MACROECONOMICS_D_L1 + MIP_FOREIGN_D_L1 +
                         ECONAPP_L1_ECONMIP_D_L1 + FPAPP_L1_FPMIP_D_L1 +
                         ECONAPP_ECONMIP_D_L1 + FPAPP_FPMIP_D_L1 + ECONAPP_L1_ECONMIP_D + FPAPP_L1_FPMIP_D +
                         APPROVE_L1 + PARTY_IN + PARTY_OUT + UNRATE_D + 
                         PARTY_IN_L1 + PARTY_OUT_L1 + UNRATE_D_L1 + 
                         DIVIDEDGOV + ELECTION + HONEYMOON +
                         ghwb + bc + gwb + bo,
                       data = approval_standard)
summary(cavari.model.differenced)

# Following our own advice, we will test to see whether any parameters can be safely omitted via LASSO
polywog.fit <- polywog(APPROVE ~ APPROVE_ECONOMY + APPROVE_FOREIGN + MIP_MACROECONOMICS_D + MIP_FOREIGN_D +
                         ECONAPP_ECONMIP_D + FPAPP_FPMIP_D + 
                         APPROVE_ECONOMY_L1 + APPROVE_FOREIGN_L1 + MIP_MACROECONOMICS_D_L1 + MIP_FOREIGN_D_L1 +
                         ECONAPP_L1_ECONMIP_D_L1 + FPAPP_L1_FPMIP_D_L1 +
                         ECONAPP_ECONMIP_D_L1 + FPAPP_FPMIP_D_L1 + ECONAPP_L1_ECONMIP_D + FPAPP_L1_FPMIP_D +
                         APPROVE_L1 + PARTY_IN + PARTY_OUT + UNRATE_D + 
                         PARTY_IN_L1 + PARTY_OUT_L1 + UNRATE_D_L1 + 
                         DIVIDEDGOV + ELECTION + HONEYMOON +
                         ghwb + bc + gwb + bo,
                       data = approval_standard, unpenalized = c("ghwb", "bc", "gwb", "bo"),
                method = "alasso", family = "gaussian", degree = 1, nfolds = 10)
summary(polywog.fit)

# Interactive structure implied from polywog. A mess of cross-lags
#             APP          MIP_D           INTAXN
# ECON        t, t-1	   no, no			 t/t-1, t-1/t
# F           t, t-1       no, no          t-1/t, t-1/t-1

#Refit with selected parameters. we will keep the general structure but pare down the other X
######### Table SM. 2; ``Restricted'' model #########
polywog.model <- lm(APPROVE ~ APPROVE_ECONOMY + APPROVE_FOREIGN + MIP_MACROECONOMICS_D + MIP_FOREIGN_D + # all t terms
                  APPROVE_ECONOMY_L1 + APPROVE_FOREIGN_L1 + MIP_MACROECONOMICS_D_L1 + MIP_FOREIGN_D_L1 + # all t-1 terms
                  ECONAPP_ECONMIP_D + FPAPP_FPMIP_D + # t/t
                  ECONAPP_L1_ECONMIP_D_L1 + FPAPP_L1_FPMIP_D_L1 + #t-1/t-1
                  ECONAPP_ECONMIP_D_L1 + ECONAPP_L1_ECONMIP_D + #t/t-1 and t-1/t 
                  FPAPP_FPMIP_D_L1 + FPAPP_L1_FPMIP_D + #t/t-1 and t-1/t
					APPROVE_L1 + PARTY_IN + PARTY_OUT_L1 +
                    ghwb + bc + gwb + bo, data = approval_standard)

summary(polywog.model)

